#include "yocto/utest/run.hpp"

#include "core/verlet.hpp"
#include "core/fcc.hpp"
#include "core/lj.hpp"
#include "core/ljc.hpp"
#include "core/constants.hpp"
#include "core/full-dyn.hpp"

#include "yocto/auto-ptr.hpp"
#include "yocto/ios/ocstream.hpp"

using namespace FastMD;

static inline void save( const char *name, bool isothermal, const Frame &F, const Real t, bool append )
{
	string filename = name;
	if( isothermal )
		filename += "-nvt.dat";
	else 
		filename += "-nve.dat";
	
	ios::ocstream fp( filename, append );
	fp("%g %g %g %g\n", t, F.V, F.K, F.V + F.K );
	filename = name;
	filename += ".xyz";
	F[ Frame::POSITIONS].save_xyz( filename, "H", vformat("t=%g", t), append );
}

static inline void save_all( const char *name, bool isothermal, const Frame &F, const Real t, bool append )
{
	string filename = name;
	if( isothermal )
		filename += "-nvt.gz";
	else 
		filename += "-nve.gz";
	F.save(filename, t, append);
}

YOCTO_UNIT_TEST_IMPL(sim)
{
	Constants::Reduced();
	const size_t m  = 4;
	const U3D    cr(m,m,m);
	const size_t np = cr.items() << 2;
	const Real   dt = 0.0005;
	auto_ptr<Step>   step( new Verlet(dt) );
	Frame        F( step->vectors(), np );
	const Real   L     = FASTMD_F(1.1) * __pow( np, __ONE/__THREE );
	F.set_box( V3D(L,L,L) );
	
	FCC::Build( F[ Frame::POSITIONS ], cr, F.box, F.ibox, 0 );	
	auto_ptr<Potential> potential( new LJ( 1, 1 ) );
	bool isothermal = false;
	for( int i=1; i < argc; ++i )
	{
		const string opt = argv[i];
		if( opt == "cutoff" )
			potential.reset( new LJCutOff(1,1,2.5) );
		if( opt == "isothermal" )
		{
			std::cerr << "ISOTHERMAL!" << std::endl;
			isothermal = true;
		}
	}
	auto_ptr<Dynamics>  dynamics(  new FullDynamics() );
	
	potential->save( "./", 0.9, 3 );
	
	Real         t     = 0;
	const Real   tmax  = 10;
	const size_t niter = size_t(ceil( tmax / dt ));
	
	//-- initialize
	for( size_t i=0; i < np; ++i ) 
	{
		const Real mass = F[ Frame::VELOCITIES ][i].data = __ONE + 4*Rand();
		F[ Frame::POSITIONS ][i].data = __ONE / mass;
	}
	F.initialize_velocities( 1.0 );
	std::cerr << "T/Frame=" << F.K / ( 1.5 * R * np ) << std::endl;
	
	dynamics->eval( F, *potential, isothermal );
	std::cerr << "Initial V=" << F.V << std::endl;
	save(     potential->name(), isothermal, F, t, false );
	save_all( potential->name(), isothermal, F, t, false );
	
	Real t_old = 0;
	for( size_t iter=1; iter <= niter; ++iter )
	{
		step->predict( F );
		dynamics->eval( F, *potential, isothermal );
		step->correct( F, isothermal);
		dynamics->eval( F, *potential, isothermal );
		t = dt * iter;
		if( t - t_old >= 0.01 )
		{
			save( potential->name(),     isothermal, F, t, true );
			save_all( potential->name(), isothermal, F, t, true );
			std::cerr << "t=" << t << " V=" << F.V/np << " K=" << F.K/np << " E=" << (F.K + F.V)/np << "        \r";
			std::cerr.flush();
			t_old = t;
		}
	}
	std::cerr << std::endl;
	
	
	
	
}
YOCTO_UNIT_TEST_DONE()
